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ABSTRACT 

The  global  behaviour  of  a  class  of  predator-prey  systems,  modelled  by  a 
pair  of  non-linear  ordinary  differential  equations,  under  constant  rate 
harvesting  and/or  stocking  of  both  species,  is  presented.  Theoretically 
possible  structures  and  transitions  are  developed  and  validated  by  computer 
simulations.  The  results  are  presented  as  transition  loci  in  the  F-G  (prey 
harvest  rate  -  predator  harvest  rate)  plane. 
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SIGNIFICANCE  AND  EXPLANATION 


This  report  gives  general  and  complete  results  on  the  coexistence  of 
predator-prey  systems  under  constant  rate  harvesting  and  stocking.  A 
theoretical  analysis  is  presented  and  a  number  of  computer  simulations  are 
included.  Certain  phenomena  not  included  in  the  analysis  are  uncovered. 

The  significance  of  the  work  is  in  showing  that  different  structures  of 
the  system  can  exist ,  that  coexistence  regions  decrease  as  harvest  rates 
increase,  and  that  contrary  to  indications  given  by  system  linearisation, 
collapse  can  occur  well  before  mathematically  critical  harvest  rates  are 
reached.  This  information  should  be  useful  in  resource  management  programs. 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


COEXISTENCE  PROPERTIES  OF  SOME  PREDATOR-PREY 
SYSTEMS  UNDER  CONSTANT  RATE  HARVESTING  AND  STOCXING* 


1  2 
F.  Brauer  and  A.  C.  Soudack 


1  .  Introduction 

In  a  sequence  of  papers  [Brauer,  Soudack  and  Jarosch  (1976);  Brauer  and  Soudack 
(1979a);  Brauer  and  Soudack  (1979b);  Brauer  and  Soudack  (1980)],  we  have  analyzed  the 
qlobal  behaviour  of  predator-prey  systems  under  constant  rate  harvesting  of  either  species 
and  under  constant  rate  stocking  (which  may  be  viewed  as  negative  harvesting)  of  either 
species  and  of  both  species  simulataneously.  In  all  of  this  previous  work,  the  focus  of 
our  attention  has  been  on  the  nature  of  the  phase  portrait  for  a  given  harvest  rate  and  on 
the  transitions  between  types  of  behaviours  as  the  harvest  rate  is  changed. 

In  this  paper  we  generalize  our  earlier  results  to  allow  two  independent  constant 
harvest  rates  (positive  or  negative)  for  the  two  species.  We  also  make  a  slight  change  in 
emphasis  and  concentrate  on  classification  of  regions  in  the  harvest-rate  plane  rather  than 
on  the  phase  portraits.  Of  course,  the  phase  portrait  analysis  can  also  be  carried  out, 
and  some  such  analysis  is  necessary  for  the  desired  classification  of  regions.  However,  we 
have  chosen  to  suppress  much  of  this  detail  in  order  to  concentrate  on  the  question  of 
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what  qualitative  behaviour  ia  to  be  expected  under  changes  of  one  or  both  harveat  rates. 
Even  though  the  class  of  models  studied  is  unrealistically  simple,  our  results  imply 
warnings  about  unexpected  dangers  which  may  be  relevant  to  the  management  of  real-life 

systems . 


x  -  X  f (x,y)  -  F 

I 

y  ■  y  g<x,y)  -  g 

as  a  model  for  the  sizes  x(t)  of  a  prey  population  and  y(t)  ^f  a  predator  population. 
Here,  f(x,y)  and  g(x,y)  are  the  respective  per  capita  growth  rates  of  the  two 
population  sizes.  As  In  our  previous  work  [Brauer,  Soudack  and  Jarosch  (1976);  Brauer  and 
Soudack  (1979a);  Brauer  and  Soudack  (1979b);  Brauer  and  Soudack  (1980)]  we  assume  that 
these  depend  only  on  the  population  sizes  at  time  t.  The  prey  species  is  harvested  at  a 
constant  time  rate  F,  while  the  predator  species  is  harvested  at  a  constant  time  rate 
G.  we  permit  either  F,  or  G,  or  both,  to  be  negative,  to  represent  stocking  rather 
than  harvesting  of  the  corresponding  species. 

The  predator-prey  nature  of  the  model  is  expressed  by  the  assumptions 

(2)  fx(x,y)  <  0,  gH(x,y)  >  0,  gy(x,y)  <  0 

for  x  >  0,  y  >  0  (subscripts  denoting  partial  derivatives.  These  assumptions  imply  that 
the  equation  g(x,y)  -  0,  representing  the  (unharvested)  predator  isocline,  defines  x  as 
a  monotone  non-decreasing  function  x  -  T(y)  for  0  <  y  <  «.  We  assume  that  this  isocline 
intersects  the  x-axis  at  (J,0);  that  is,  that  J  •  HO),  or 

(3)  g< J,0)  -  0  . 

As  we  have  pointed  out  previously  (Brauer  and  Soudack  (1979a)),  in  many  of  the  commonly- 
used  models  for  function  g  is  independent  of  y,  corresponding  to  the  situation  in  which 
the  predators  do  not  interfere  with  one  another  in  their  search  for  prey.  In  this  case  the 
predator  Isocline  g(x,y)  “0  is  the  vertical  line  x  =  J. 

The  hypothesis  (2)  also  implies  that  the  equation  f(x,y)  “  0,  representing  the 
(unharvested)  prey  isocline,  defines  y  as  a  single-valued  function  y  ”  4(x)  which  we 
assume  non-negative  on  an  interval  0<o<x<K<*  with  4(K)  »  0,  or 

(4)  f(K,0)  -  0  . 


2.  General  Theory 

We  consider  the  system 

(1) 


-3- 


It  la  necessary  to  distinguish  three  possibilities,  as  follower 
(i)  a  >  0,  corresponding  to  f(0,0)  <  0. 

(ii)  a  -  0,  corresponding  to  f(0,0)  >  0,  and  there  exists  L  <  "  such  that 

(5)  f<0,L>  -  0 

(ili)  a  -  0,  L  -  -  . 

Biologically,  a  >  0  is  the  case  in  which  the  prey  population  is  unable  to  develop  if  it 

gets  too  small,  even  in  the  absence  of  predators.  If  a  -  0,  the  prey  population  can 
establish  itself  front  a  small  initial  population;  L  is  the  maximum  predator  density  for 
which  the  prey  population  can  establish  itself.  If  a  »  0,  L  «  »,  the  prey  population  can 
establish  itself  for  any  predator  population.  Finally,  we  assume 

(6)  a  <  J  <  X  . 

The  cases  J  >  K  and  J  <  a  may  be  analyzed  by  the  same  methods  and  are  essentially 
trivial.  A  discussion  of  the  biological  significance  of  the  numbers  a,  J,  K,  L  may  be 
found  in  [Brauer  and  Soudack  (1979a)]. 

*  *  a 

An  equilibrium  P(x,y)  of  the  system  (1)  is  an  intersection  of  the  prey  isocline 
xf (x,y)  -  F  *  0  and  the  predator  isocline  yg(x,y)  -  G  «  0.  To  study  the  (local) 
stability  of  an  equilibrium,  we  linearize  about  the  equilibrium,  forming  the  matrix 

*  f  (x,y)  +  f(x,y)  x  f  (x,y) 

A(p)  -  x  a  y 

Y  9„<x,y>  y  g  (x,y)  +  g(*,y> 

*  y 

and  then  determine  the  eigenvalues  of  A(P).  In  particular,  if  det  A(P)  <  o,  then  the 
eigenvalues  have  opposite  sign  and  P  is  a  saddle  point.  If  det  A(P)  >  0,  then  the  real 
parts  of  the  eigenvalues  have  the  same  sign  and  P  is  a  node  or  spiral  point  which  is 
asymptotically  stable  if  tr  A(P)  <  0  and  unstable  if  tr  A(P)  >  0  (corresponding  to 
eigenvalues  with  negative  real  part  and  positive  real  part  respectively) . 

If  GM,  the  predator  isocline  y  g(x,y)  -  G  is  a  curve  which  approaches  the 
curve  g(x,y)  »  0  asymptotically.  The  portion  of  this  curve  in  the  first  quadrant  lies  to 
the  right  of  g(x,y)  -  0  if  G  >  0  and  to  the  left  of  g(x,y)  -0  if  G  <  0. 


To  describe  the  prey  isocline  x  f(x,y)  -  F  for  F  ^  0,  we  define 

*  F 

f  <x,y)  «  f(x,y)  -  -  • 

Then  the  prey  isocline  is  the  curve  f*(x,y)  »  0;  we  shall  regard  f  (x,y)  as  a  modified 

per  capita  growth  rate  for  each  fixed  F.  it  is  easy  to  see  that  for  every  F  <  0,  f 

satisfies  the  sane  hypotheses  as  f  and  is  of  the  type  a  «  0,  L  *  independent  of  the 

type  for  F  -  0.  For  F  >  0,  there  are  two  critical  values  of  F.  For  sufficiently 

small  F  >  0,  f*  satisfies  the  same  hypotheses  as  f  and  is  of  the  type 

o  >  0,  with  a  «  a(F)  <  J  and 

K  -  K(F)  >  J.  There  exists  Fc  >  0  such  that  either 

(i)  o(F  )  <  J,  K(F  )  -  J 
c  c 

or 

(ii)  o(F  )  -  J,  K(F  )  >  J  . 
c  c 

Further,  there  exists  F*  >  Fc  for  which 

a(F*>  -  K(F*) 

[Brauer  and  Soudack  (1979b)].  Then  f*  satisfies  the  same  hypotheses  as  f  and  is  of  the 
type  o  >  0  for  0  <  F  <  F*. 

This  suggests  that  models  of  the  type  a  >  0  with  F  -  0  may  be  the  result  of  some 
harvesting  of  prey  by  biological  mechanisms  outside  the  system,  while  models  of  the  type 
a  «  0,  L  «  •  may  represent  some  external  stocking  of  prey.  Models  with  a  «  0,  L  <  “ 
which  are  the  type  most  frequently  studied,  represent  a  critical  balance  between  external 
harvesting  and  stocking. 

For  each  suitably  chosen  fixed  F,  we  may  now  view  the  system  (1)  as  a  pure  predator 
harvesting  or  stocking  model 

•  » 
x  ■  x  f  (x,y ) 

(7) 

y  -  y  g(x,y)  -  g  , 

with  the  prey  harvest  or  stocking  built  into  the  prey  growth  rate.  Thus  as  G  is  varied 
for  any  fixed  F,  the  possible  transitions  are  of  the  same  types  as  for  pure  predator 
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harvesting  (Brauer  and  Soudack  ( 1979a) ) •  The  first  question  we  wish  to  examine  is  the 


existence  of  an  equilibrium  of  the  system  (1),  or  (7),  in  the  interior  of  the  first 

quadrant.  It  is  known  [Brauer  and  Soudack  (1979a)]  that  there  exists  G+  «  g+(f)  >  0  such 

that  the  system  (1)  has  an  equilibrium  PK  •  Pm(F,G),  not  a  saddle  point,  in  the  interior 
of  the  first  quadrant  of  the  x-y  plane  if  0  <  G  <  G+,  and  has  no  equilibrium  in  the 

interior  of  the  first  quandrant  if  G  >  G+.  This  holds  for  all  F  <  0  and  for 

»  * 

sufficiently  small  positive  F.  It  is  easy  to  see  graphically  that  G  (F)  must  be  a 
monotone  decreasing  function  of  F  for  -“>  <  F  <  ". 

If  G  <  0  we  must  distinguish  between  the  cases  F  >  0  and  F  <  0  because  of  the 
difference  between  the  cases  a  >  0  and  a  “  0,  L  «  °»,  [Brauer  and  Soudack  (1990)].  If 
F  >  0,  corresponding  to  a  >  0,  there  exists  G"  ■  G”(F)  <  0  such  that  the  system  (1) 

has  an  equilibrium  P^  *  P^tFjG)  (not  a  saddle  point)  in  the  interior  of  the  first 

quadrant  if  0  >  G  >  G  ,  and  has  no  equilibrium  in  the  interior  of  the  first  quadrant  if 
G  <  G”.  It  is  easy  to  see  that  G“(F)  is  a  monotone  incresing  function  of  F.  If  F  < 

0,  corresponding  to  a*  0,  !■•*,  the  system  (1)  haB  an  equilibrium  P^  in  the  interior 
of  the  first  quadrant  for  all  G  <  0. 

The  above  analysis  is  valid  for  all  F  <  0,  but  for  positive  F  it  is  obviously 

valid  only  for  0  <  F  <  F  .  However,  if  a(F  )  <  J,  K(F  )  «  J  and  G  <  0,  it  is  easy  to 

c  c  c 

see  that  it  is  in  fact  valid  for  0  <  F  <  F  .  Similarly,  if  a(Fc>  *  J,  K(FC)  >  j,  and  G 

* 

>  0,  it  is  valid  for  0  <  F  <  F  .  We  may  collect  this  information  on  the  values  of  the 
harvest  rates  F  and  G  for  which  there  is  an  equilibrium  P^  «  P^FjG)  in  the  interior 
of  the  first  quadrant  of  the  x-y  plane.  In  Figures  1  and  2,  the  interior  of  the  shaded 
region  indicates  the  set  R  of  such  values  in  the  F-G  plane.  Such  figures  have  been 
given  before  for  species  in  competition  [Yodzis  (1976),  Reading  (unpublished),  Griffel 


(1979)] 


As  we  have  shown  in  our  earlier  work  on  harvesting  of  one  species  [Brauer  and  Soudack 
(1979a,  1979b)],  the  existence  of  an  equilibrium  in  the  interior  of  the  first  quadrant  of 
the  x-y  plane  does  not  guarantee  the  survival  of  both  species.  Resolution  of  this 
question  depends  not  only  on  the  existence  of  an  equilibrium  but  also  on  the  structure  of 
separatrices  at  the  one  or  more  saddle  points  in  the  interior  of  the  first  quadrant.  We 
have  given  a  classification  in  the  two  one-species  harvesting  problems  which  can  be 
combined  into  the  following  classification  for  any  model  of  the  form  (1)  with  (F,G)  in 
the  interior  of  R,  so  that  there  is  an  equilibrium  P^tF, G)  In  the  interior  of  the  first 
quadrant  of  the  x-y  plane.  There  is  then  at  least  one  saddle  point  in  the  first  quadrant, 
possibly  on  one  of  the  axes. 

Case  1.  There  is  an  orbit  running  from  a  saddle  point  as  t  ♦  -»  to  or  a  limit 

cycle  around  p  as  t  ♦  +*•. 

Case  2.  There  is  an  orbit  running  from  a  saddle  point  as  t  *  -«•  to  a  saddle  point 
as  t  ♦  +•  (homoclinic-type  orbit).  The  two  saddle  points  may  be  the  same. 

Case  3.  There  is  an  orbit  running  from  PM  or  a  limit  cycle  around  as  t  +  -«• 

to  a  saddle  point  as  t  *  +". 

For  each  case  there  are  two  alternatives  which  we  index  as  a  or  b  according  as  the 
equilibrium  P^  is  (locally)  asymptotically  stable  or  unstable  respectively.  In  case  la, 
P^  is  asymptotically  stable  and  there  is  a  region  of  asymptotic  stability  for  PM  -  the 
set  of  initial  values  for  which  the  solution  tends  to  Pm  as  t  *  ■>,  which  can  be 
described  in  terms  of  the  separatrices  at  saddle  points.  In  case  1b,  P^  is  unstable  but 
there  is  an  asymptotically  stable  limit  cycle  with  a  domain  of  asymptotic  stability  which 
again  can  be  described  in  terms  of  separatrices.  Case  2  may  be  viewed  as  a  transition 
between  case  1  and  case  3.  In  case  3a,  P^  is  asymptotically  stable,  but  the  domain  of 
asymptotic  stability  consists  only  of  the  interior  of  an  unstable  periodic  orbit  around 
P„.  In  case  3b,  P,,  is  unstable,  and  every  orbit  goes  to  an  axis  in  finite  time, 

corresponding  to  extinction  of  one  of  the  species.  Thus  in  practical  terms,  the  two 
species  can  coexist  only  if  the  system  is  in  case  1,  even  though  in  case  3a  there  is  a 


"small”  set  of  Initial  values  for  which  there  is  coexistence,  and  in  case  2  there  may  be  a 
coexistence  region  which  is  extremely  susceptible  to  collapse  under  small  perturbations  of 
harvest  rates. 

This  suggests  the  importance  of  dividing  the  region  R  in  the  F-G  plane  into 
subregions  corresponding  to  the  various  cases  in  order  to  determine  the  set  of  values 
(F,G)  for  which  the  system  can  continue  to  function  with  both  species  co-existing.  For 
any  given  (F,G)  e  R,  we  may  calculate  Pm(F,G),  and  then  tr  MP^fF.G)}.  The  set  of 
points  (F,G)  e  R  for  which  tr  i{Pai(F,G) }  •  0  is  a  curve  a  in  R  corresponding  to  the 
transition  between  the  alternatives  a  and  b  (local  asymptotic  stability  and  instability 
respectively) .  There  may  be  another  curve  h  in  R  describing  the  set  of  values 
(F,G)  e  R  for  which  there  is  a  homoclinc  type  orbit  and  the  system  is  in  case  2.  While 
the  curve  o  may  be  drawn  approximately  by  calculation  of  tr  A{Pa|(F,G)},  the  curve  h 
can  be  approximated  only  by  computer  simulation  of  the  orbits  of  the  system  ( 1 )  and 
classification  of  cases  [Brauer  and  Soudack  (1979a),  (1979b)]. 

By  examining  the  prey  and  predator  isoclines  it  is  not  difficult  to  see  that  the 
qualitative  structure  for  the  system  (1)  with  F  0  is  the  same  as  the  qualitative 
structure  for  the  system  (1)  with  F  ■  0,  except  that  there  are  differences  when  G  «  0 
as  shown  in  (Brauer  and  Soudack  (1979b)].  By  the  same  methods  as  those  used  in  the  case  F 
■  0,  we  may  establish  the  following  results. 

Theorem  1 :  There  is  a  neighbourhood  of  the  origin  in  the  F-G  plane  in  R  for  which 
the  system  (1)  is  in  case  la  or  in  case  1b. 

Theorem  2;  If  (F,G)  is  sufficiently  close  to  the  boundary  of  the  region  R  in 
the  F-G  plane  for  which  the  system  (1)  has  an  equilibrium,  and  if  tr  A{Pw(F,G)}  f  0, 
then  the  system  (1)  is  in  case  la  or  in  case  3b. 

It  follows  from  Theorem  2  that  if  the  curve  h  goes  to  the  boundary  of  R,  it  must 
intersect  the  curve  a  there,  since  h  can  meet  the  boundary  of  R  only  in  a  point  where 
tr  ) }  •  0.  We  have  also  shown  [Brauer  and  Soudack  (1980)]  that  in  the  third 

quadrant  of  the  F^S  plane,  corresponding  to  stocking  of  both  species,  the  system  is  in 
case  la  or  case  1b.  Thus  the  curve  h  can  not  enter  the  third  quadrant. 


A  final  general  remark  ie  that  aa  we  have  ahown  in  [Brauer  and  soudack  (1979b)},  if 
the  boundary  G  -  G+(y)  of  R  meets  the  F-axia  at  <Fe.O)  and  the  bourdary  G  •  G"(F) 
meets  the  F-axia  at  (F* ,0 )  (Figure  1),  then  the  system  is  in  ease  la  at  t  Fc , 0 ) ,  while 
if  G  “  G+(F)  meets  the  F-axis  at  (F*,0)  and  G  ■  G"(F)  meets  the  F-axle  at  (Fc»0), 
(Figure  2)<  then  the  system  is  in  case  3b  at  (Fc,0). 

In  the  next  section  we  shall  indicate  by  considering  a  class  of  examples  how  the 
classifications  for  pure  predator  harvesting  and  pure  prey  harvesting  can  give  information 
about  the  structure  of  the  region  R  and  the  claesification  for  two  -  species  harvesting. 


3.  A  Class  of  Examples 


In  our  previous  work  on  harvesting  and  stocking  [Brauer,  Soudack  and  Jarosch  (1976); 
Brauer  and  Soudack  (1979a),  (1979b),  (1980)],  we  have  used  the  model 

f(x,y)  -  r(l  -  f)  -  Jfj 

(8) 

,  >  ,  x  J  .  sA(x-J) 

g(x,y)  -  8<X+A  ~  J+A>  "  (j+A)(x+A) 

[Holling  (1965)]  as  a  source  of  examples.  For  this  model,  K  and  J  are  the  K  and  J 
of  the  general  theory  in  Section  2,  while  L  =  rA.  As  we  have  shown  previously. 
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It  is  easy  to  calculate  that 


tr  A(x,y)  -  r(1  -  — )  -  ■ 

In  particular,  if  F  «  G  »  0,  then  xK  »  J,  y( 

tr  £1x^(0, 0) ,  y^tOjO)] 

For  F  «  Fc  and  G  -  0,  we  have  x<>  »  J,  yw  • 

tr  AlxJF^O),  yjlc 

For  F  »  0  and  G  -  Gc,  we  have  xm  «  21*.  y( 

tr  4[x  (0,G  ),  y  (0,G  )]  - 
«  c  00  c 


Ay 


(x+A) 


sA(  x-J) 
(J+A) (x+A) 


-  |  (J+A)  (K-J)  ,  and 


rJ 

K( J+A) 


(K-A-2J) 


0,  and 


,0)]  -  |  (K-2J)  . 


_  r  ( J+K+2A) (K-J) 

K  -  4 - ’  and 


8AK(K-J)-r(A+J)2(J+K) 
K( J+A) (2A+J+K) 


The  stability  curve  o  in  the  (F,G)  plane  corresponds  to  tr  A  »  0.  Thus  for  a 
given  pair  (F,G),  we  calculate  P^tF.G),  and  hence  tr  &[PC<(F,G)].  By  varying  (F,G) 
we  may  plot  the  curve  a.  The  curve  h  describing  the  pairs  (F,G)  for  which  there  is  a 
homoclinic  orbit  can  n.-.t  be  sketched  so  easily.  For  this  we  need  information  about  the 
dynamics  of  the  system,  and  must  compute  orbits. 

As  we  have  observed  in  the  general  theory,  there  are  situations  in  which  we  may  have 
either  predator  extinction  or  prey  extinction,  depending  on  the  initial  Btate.  We  remark 
that  it  is  easy  to  separate  these  two  possibilities  by  computing  the  orbit  from  the  orioin 
in  the  x-y  plane  backwards  in  time.  This  orbit  serves  as  a  separatrix  between  predator 
extinction  and  prey  extinction.  (In  practice,  the  computation  may  be  carried  out  more 
efficiently  if  the  starting  point  is  taken  near  the  origin  rather  than  at  the  origin. 

We  now  give  several  examples  using  this  model  to  indicate  the  range  of  possibilities 
and  the  procedure  for  analyzing  a  model.  For  a  given  choice  of  the  parameters  of  the  model 
and  a  given  pair  of  harvest  rates,  orbits  can  be  approximated  by  computer  simulation  just 
as  in  our  previous  work  [Brauer,  Soudack  and  Jarosch  (1976),  Brauer  and  Soudack  (1979a, 
1979b,  1980)],  in  this  paper,  we  emphasize  the  shape  and  structure  of  the  region  R  in 
the  F-G  plane  for  which  there  is  an  equilibrium  in  the  interior  of  the  first  quadrant  of 
the  x-y  plane,  the  curve  o  in  R  corresponding  to  the  transition  between  local 
asymptotic  stability  and  instability  of  this  equilibrium,  and  the  curve  h  in  R 
corresponding  to  values  (F,G)  for  which  there  is  a  homoclinic  type  orbit.  Thus  we  show 
orbits  for  only  one  of  the  examples,  even  though  it  is  necessary  to  compute  orbits  in  order 
to  locate  the  curve  h  for  every  example.  The  curve  a  and  the  boundary  of  the  region 
R  can  be  approximated  numerically  by  calculation  of  equilibria  and  the  trace  of  the  matrix 
A,  without  requiring  calculations  of  orbits.  The  computations  reported  below  were 
carried  out  on  the  University  of  Wisconsin  UNIVAC  1110  and  the  University  of  British 
Columbia  Amdahl  470. 

Our  general  procedure  in  each  example  has  been  to  examine  the  case  transitions  as  one 
of  the  harvest  rates  F,G  is  varied  while  the  other  is  held  at  zero.  This  gives  the 
"corners"  of  the  region  R.  Then  the  boundary  of  R  can  be  filled  out  by  varying  G 


with  F  fixed.  To  locate  the  curvea  a  and  h,  we  proceed  juat  aa  we  did  in  our 
previoua  work  dealing  with  harveeting  and  stocking  of  a  single  species.  For  each  example 
we  make  some  general  observations,  not  all  of  which  are  shown  in  the  figures  reproduced 
here. 

Example  li  r  •  1,  s  ■  1,  A  ■  10,  J  »  20,  K  ■  40. 

For  this  set  of  parameters,  Fc  *  F*  -  10,  Ge  -  0.83,  from  (9).  The  system  is  in  case  la 
for  all  (F,G)  e  A  (Figure  3)  and  is  strongly  stable  in  the  sense  that  orbits  approach  the 
equilibrium  Pv(F,G)  rapidly.  The  orbits  are  similar  to  those  in  the  case  of  pure 
predator  harvesting  in  that  the  stability  region  determined  by  the  separatrlces  at  the 
saddle  point  becomes  smaller  as  G  is  increased,  for  each  F.  Outside  the  region  of 
asymptotic  stability,  there  is  a  region  of  prey  extinction  and  a  region  of  predator 
extinction  which  may  be  separated  numerically  by  integration  of  the  system  backwards  in 
time  from  (c,c). 
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Fig.  3 


Example  2;  r  •  1,  s  "  4,  A  •  10,  J  -  20,  X  «  60. 

Here,  Gc  -  8.88,  rc  -  13.33.  F*  -  15,  from  (9).  In  ooth  F  and  G  separately,  the  case 
transitions  are  1b  ♦  2b  ♦  3b  (Figure  4).  Note  that  since  the  equilibrium  is  unstable 
at  <Fc,o)  the  region  R  extends  to  <F*,0)  in  the  first  gradrant.  As  in  Example  1,  the 
region  of  asymptotic  stability  becomes  smaller  as  G  increases,  and  there  is  a  coexistence 
region  only  for  a  small  part  of  R. 
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Fig.  4 


Example  3:  r  ■  1,  s  ■  7,  A  »  10,  J  •  20,  K  »  40. 

Here,  G  -  5.83,  F_  “  P*  »  10.  The  caee  transition  is  la  +  2a  ♦  3a  *  3b  in  G  for  P  * 

C  c 

* 

0,  and  the  system  is  in  cese  la  for  G  ■  0,  0<F<F.  The  3a  region  is  very  small, 
suggesting  that  in  practical  terms  the  transition  is  la  *  3b.  The  system  is  weakly  stable 


inn  the  sense  that  orbits  approach  the  equilibrium  tm  very  slowly.  For  this  example,  we 
have  also  indicated  the  region  R  in  all  four  quadrants  of  the  F-G  plane  (Figure  5). 
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Example  4:  r  ■  1,  s  »  5,  A  -  10,  J  -  20,  K  -  45. 

Here,  Gc  “  5.787,  Fc  «  11.lT,  F*  -  11.25.  The  caee  transltlone  are  la  +  1b  *  2b  ♦  3b 
in  F  and  la  ♦  1b  ♦  la  in  G.  Figure  6  indicates  that  the  curves  o  and  h  meet  on 
the  boundary  of  the  region  R.  There  are  two  components  of  the  curve  o  in  the  first 
quadrant;  in  Figure  6,  we  see  how  the  curve  o  goes  into  the  second  quadrant  and  then 
returns  to  the  first  quadrant,  so  that  the  two  components  are  in  the  first  quadrant  are  not 
really  separate.  The  reader  should  note  from  Figure  6  that  for  fixed  F  •  1  and 
increasing  G,  the  case  transition  is  la  ♦  1b  ♦  2b  ♦  3b  ♦  2b  +  1b  ♦  1a»  this  is  much 
more  complicated  than  for  F  -  0,  but  can  still  be  read  from  the  figure. 
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Example  5;  r  -  2,  a  ■  1,  A  ■  10,  J  -  20,  K  »  60. 

Here,  G  “  4.44,  F  ■  26.6,  F*  “  30.  The  case  transitions  are  1b  ♦  la  in  G  and 
c  c 

1b  +  2b  ♦  3b  in  F.  Again,  as  may  be  seen,  several  transitions  are  possible  when  one 

of  F  and  G  is  held  fixed  and  the  other  is  varied,  more  than  in  the  special  cases  F  » 

0  and  G  •  0.  US  have  included  sane  phase  portraits  for  this  example.  He  see  (Figure  7) 
that  for  F  ■  3,  G  «  0.5  the  system  is  in  case  3b.  For  F  -  3,  G  «  1  (Figure  8)  the 
system  has  just  shifted  to  case  Ibj  the  limit  cycle  is  large  with  a  further  increase  in 

G  to  F  •  3,  G  ■  2  (Figure  9)  the  system  is  still  in  case  1b  but  the  limit  cycle  is 

smaller.  For  F  «  3,  G  •  3  (Figure  10)  the  system  is  in  case  la.  Note  that  in  Figures 
11,  12,  12,  the  sepatatrix  between  predator  extinction  and  prey  extinction  is  roughly  the 


same  but  the  region  of  coexistence  shrinks  as  G  increases  even  though  P^  is  stabilizing 
(locally).  For  larger  F,  F  ”  13,  G  «  1.6  (Figure  11)  and  F  ■  13,  G  “  2  (Figure  12)  we 
have  cases  1b  and  la  respectively.  The  predator  extinction  -  prey  extinction  separatrix  is 
lower,  giving  a  larger  region  of  prey  extinction  and  smaller  regions  of  coexistence  and 
predator  extinction.  In  Figure  13,  we  have  shown  all  four  quadrants  of  the  F-G  plane, 
with  the  curves  o  and  h  meeting  at  two  points  of  the  boundary  of  R.  The  1b  region  in 
the  fourth  quadrant  is  unstable  in  practical  terms  because  the  limit  cycles  come  very  close 
to  the  coordinate  axes.  There  is  also  a  small  la  region  in  the  fourth  quadrant,  but  this 
also  is  not  practically  stable  because  P  is  very  close  to  the  y-axis. 
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C-COEXISTENCE  REGION 


Fig.  13 


We  have  also  examined  some  variations  on  Example  5,  holding  r  «  2, 
a  “  1,  A  ■  10,  J  «  20  and  changing  K.  For  example,  with  JC  <*  50,  the  curve  a  passes 
through  the  origin  (an  immediate  consequence  of  the  fact  that 
K  ■  2 J  +  A) .  Qualitatively,  the  picture  (Figure  14)  is  much  like  that  for 
K  «  60  (Figure  13),  except  that  the  curves  a  and  h  have  moved  down.  Similarly,  for  K 
-  70  the  picture  is  similar  except  that  the  curves  a  and  h  move  upwards. 
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4.  Conclusions. 

In  this  paper,  we  have  generalized  and  completed  the  study  of  constant-rate  harvesting 
and  stocking  in  a  class  of  idealized  predator-prey  systems.  In  general,  the  existence  of  a 
saddle  point  in  the  first  quadrant  of  the  phase  plane  causes  the  trajectories  to  be  similar 
to  those  for  predator  harvesting  alone  rather  than  to  prey  harvesting  alone  (see  References 
2  and  3).  There  is  however,  one  more  transition  trajectory  separating  the  regions  of 
predator  extinction  and  prey  extinction,  as  would  be  surmised. 

Thus,  we  confined  our  attention  to  the  F-G  plane,  and  utilizing  the  fact  of 
continuity  across  the  F  and  G  axes,  we  established  a  region  R,  for  which  exists 

and  loci  a  and  h  which  divide  R  into  regions  of  stability,  local  instability,  and  no 
coexistence  of  species  at  all.  The  behaviour  of  the  system,  for  all  F  and  G,  including 
all  combinations  of  harvesting  and  stocking,  is  described  by  the  loci  in  this  plane. 

One  notable  difference  between  two-species  and  one-species  harvesting  is  that  we  no 
longer  must  stay  in  Case  3b  once  we  get  there.  Increasing  harvest  rates  might  again 
establish  coexistence  regions.  This  phenomenon  does  not  exist  if  either  species  is 
harvested  separately. 

The  case  of  two-species  stocking  is  not  examined  in  detail  here,  since  the  results  are 
given  in  [Brauer  and  Soudack,  1980].  It  is  enough  to  notice  that  the  a  locus  never 
enters  the  3r<^  quadrant  of  the  F-G  plane  and  thus  there  is  a  region  of  coexistence  for 
all  stocking  rates. 

Examples  where  one  specie  is  stocked  and  the  other  harvested,  corresponding  to  the  2n<1 
and  4th  quadrants  of  the  F-G  plane,  were  considered,  but  the  results  are  not  included 
because  the  phase  portraits  are  merely  combinations  of  the  cases  for  single  specie 
harvesting  and  stocking.  The  relevant  information  is  included  in  the  description  of  the 
F-G  plane. 

Admittedly,  there  is  some  question  regarding  the  idealized  models  considered  in  this 
study.  However,  studies  of  these  idealized  models  have  shown  that  the  results  obtained  by 
linearization  around  equilibrium  points  are  often  incorrect  and  that  such  systems  tend  to 
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be  much  mere  fragile  than  expected,  aa  far  aa  coexistence  is  concerned.  The  next  step  is 
to  consider  more  realistic  harvesting  strategies  and  to  expand  the  models  to  include  age 
and/or  sice  structure. 

The  authors  wish  to  thank  Mr.  A1  MacKencie  of  the  Department  of  Electrical 
Engineering,  University  of  British  Columbia,  for  preparing  the  figures  in  this  paper. 
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